Resistome_visualization
Florentin Constancias
January 07, 2022
Load required packages
library(tidyverse)
library(phyloseq)
library(speedyseq)
library(ggrepel)
library(here)
library(microViz)
library(RColorBrewer)
library(vegan)
library(randomcoloR)
options(getClass.msg=FALSE) # https://github.com/epurdom/clusterExperiment/issues/66
#this fixes an error message that pops up because the class 'Annotated' is defined in two different packagesSource required functions
rm(list = ls())
'%!in%' <- function(x,y)!('%in%'(x,y))
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_taxa_tests.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_normalisation.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_alpha.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_beta.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_heatmap.R")Import data and define variables:
physeq = readRDS("~/Projects/ETH/Alessia/mobilome/ps_combined.rds")
out_pptx = "~/Projects/ETH/Alessia/mobilome/graph.pptx"
physeq %>%
tax_table() %>%
data.frame() %>%
mutate(Resistance_Mechanism_multi = ifelse(grepl(";", Resistance_Mechanism), "multi-mech", Resistance_Mechanism)) %>%
as.matrix() -> tax_table(physeq)
physeq %>%
tax_table() %>%
data.frame() -> tax_mapping
# physeq %>%
# tax_table() %>%
# data.frame() %>%
# mutate_if(is.character, as.factor) %>%
# as.matrix() -> tax_table(physeq)
#
# physeq %>%
# tax_table() %>%
# data.frame() %>%
# mutate_all(funs(str_replace(., "[^[:alnum:]]", "_"))) %>%
# mutate_all(funs(str_replace(., "[[:punct:]]", "_"))) %>%
# as.matrix() -> tax_table(physeq)
physeq %>%
otu_table() %>%
data.frame() %>%
replace(is.na(.), 0) %>%
# mutate_all(funs(str_replace(., "[^[:alnum:]]", ""))) %>%
mutate_all(as.numeric) %>%
as.matrix() %>%
otu_table(taxa_are_rows = TRUE) -> otu_table(physeq)
# df[,'values'] <- as.numeric(df[,'values'])
# otu_table(physeq) <- tmpall.equal(taxa_names(physeq %>% tax_table()), taxa_names(physeq %>% otu_table()))## [1] TRUEall.equal(taxa_names(physeq), taxa_names(physeq %>% otu_table()))## [1] TRUEphyseq_tmp <- physeq
tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Resistance_Mechanism_multi", "Best_Hit_ARO")]
# tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Best_Hit_ARO", "Drug_Class")]
physeq_tmp %>%
tax_glom(., taxrank = "Best_Hit_ARO") %>%
tax_glom(., taxrank = "Resistance_Mechanism_multi") -> physeq_res_multi
taxa_names(physeq_res_multi) <- tax_table(physeq_res_multi)[,"Resistance_Mechanism_multi"]
physeq_res_multi## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 6 taxa and 92 samples ]:
## sample_data() Sample Data: [ 92 samples by 60 sample variables ]:
## tax_table() Taxonomy Table: [ 6 taxa by 2 taxonomic ranks ]:
## taxa are rowsphyseq_tmp <- physeq
tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Drug_Class_multi", "Best_Hit_ARO")]
# tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Best_Hit_ARO", "Drug_Class")]
physeq_tmp %>%
tax_glom(., taxrank = "Best_Hit_ARO") %>%
tax_glom(., taxrank = "Drug_Class_multi") -> physeq_drug_multi
taxa_names(physeq_drug_multi) <- tax_table(physeq_drug_multi)[,"Drug_Class_multi"]
physeq_drug_multi## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 19 taxa and 92 samples ]:
## sample_data() Sample Data: [ 92 samples by 60 sample variables ]:
## tax_table() Taxonomy Table: [ 19 taxa by 2 taxonomic ranks ]:
## taxa are rowsphyseq_tmp <- physeq
tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Drug_Class", "Best_Hit_ARO")]
# tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Best_Hit_ARO", "Drug_Class")]
physeq_tmp %>%
tax_glom(., taxrank = "Best_Hit_ARO") %>%
tax_glom(., taxrank = "Drug_Class") -> physeq_drug
taxa_names(physeq_drug) <- tax_table(physeq_drug)[,"Drug_Class"]
physeq_drug## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 43 taxa and 92 samples ]:
## sample_data() Sample Data: [ 92 samples by 60 sample variables ]:
## tax_table() Taxonomy Table: [ 43 taxa by 2 taxonomic ranks ]:
## taxa are rowsphyseq_tmp <- physeq
tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Best_Hit_ARO", "Drug_Class")]
physeq_tmp %>%
tax_glom(., taxrank = "Best_Hit_ARO") -> physeq_AMRgn
taxa_names(physeq_AMRgn) <- tax_table(physeq_AMRgn)[,"Best_Hit_ARO"]
physeq_AMRgn## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 142 taxa and 92 samples ]:
## sample_data() Sample Data: [ 92 samples by 60 sample variables ]:
## tax_table() Taxonomy Table: [ 142 taxa by 2 taxonomic ranks ]:
## taxa are rowsDefine colors for Drug classes:
set.seed(123)
treat_col <- ggpubr::get_palette(palette = "lancet",
k = get_variable(physeq_drug_multi, "Treatment") %>%
levels() %>%
length())
names(treat_col) <- get_variable(physeq_drug_multi, "Treatment") %>%
levels()Define colors for Treatment - no matter concentration / fermentation / reactors:
myCol <- viridis::viridis(n = get_variable(physeq_drug_multi, "Treatment") %>%
levels() %>%
length() + 20)
myCol## [1] "#440154FF" "#470F62FF" "#481C6EFF" "#482878FF" "#463480FF" "#423F85FF"
## [7] "#3E4A89FF" "#3A548CFF" "#355F8DFF" "#31688EFF" "#2D718EFF" "#297A8EFF"
## [13] "#26828EFF" "#228C8DFF" "#1F948CFF" "#1F9E89FF" "#21A685FF" "#29AF7FFF"
## [19] "#35B779FF" "#45BF70FF" "#58C765FF" "#6DCD59FF" "#83D44CFF" "#9BD93CFF"
## [25] "#B4DE2CFF" "#CDE11DFF" "#E6E419FF" "#FDE725FF"scales::show_col(myCol,
cex_label = 0.5)
myCol <- ggpubr::get_palette(palette = "npg",
k = get_variable(physeq_drug_multi, "Treatment") %>%
levels() %>%
length() + 3)
myCol## [1] "#E64B35" "#5CAFC5" "#0FA596" "#2A6A87" "#A97E82" "#BB9599" "#89AAB9"
## [8] "#A79287" "#C9130E" "#83664E" "#B09C85"scales::show_col(myCol,
cex_label = 0.5)
myCol <- ggpubr::get_palette(palette = "jama",
k = get_variable(physeq_drug_multi, "Treatment") %>%
levels() %>%
length() + 6)
myCol## [1] "#374E55" "#846C4D" "#D28A45" "#89957B" "#229EBE" "#3685A8" "#885B66"
## [8] "#A45F57" "#8A8F7D" "#76A397" "#6F8198" "#6B6695" "#756F80" "#80796B"scales::show_col(myCol,
cex_label = 0.5)
myCol <- ggpubr::get_palette(palette = "lancet",
k = get_variable(physeq_drug_multi, "Treatment") %>%
levels() %>%
length() + 5)
myCol## [1] "#00468B" "#9E172E" "#B43C15" "#42B540" "#16A28D" "#3085AD" "#925E9F"
## [8] "#D99395" "#E2746E" "#AD002A" "#AD7987" "#7C8181" "#1B1919"scales::show_col(myCol,
cex_label = 0.5)
treat_col = get_variable(physeq_drug_multi, "Treatment") %>%
levels()
# treat_col <- c("#1B1919", "#7C8181", "#00468B", "#2589AE","#029AAF", "#463161", "#7C66A2", "#ad99cf")
treat_col <- c("#1B1919", "#7C8181", "#00468B", "#3085AD", "#42B540", "#9E172E", "#B43C15", "#AD7987")
names(treat_col) <- get_variable(physeq_drug_multi, "Treatment") %>%
levels()
scales::show_col(treat_col,
cex_label = 0.5)
set.seed(123)
col <- distinctColorPalette(ntaxa(physeq_drug_multi))
names(col) <- taxa_names(physeq_drug_multi)physeq_AMRgn %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
left_join(tax_mapping,
by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
suffix = c("_x", "")) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
separate_rows(., Drug_Class,sep = '; ',convert = TRUE) %>%
ggplot(mapping = aes(x = as.factor(Day_of_Treatment),
y = Best_Hit_ARO ,
fill = logTPM,
)) +
facet_grid(Drug_Class ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free") + #, labeller = labeller(Drug_Class_multi = label_wrap_gen(width = 15))) +
# facet_grid(. ~ Reactor_Treatment, scales = "free", space = "free", drop = TRUE,labeller = labeller(Drug.Class = label_wrap_gen(width = 15))) +
#we adjust the vertical labels such that it is printed correctly.
geom_tile() +
theme_light() +
# scale_fill_viridis_c(name = "TPM",
# na.value = "transparent", trans = scales::pseudo_log_trans(sigma = 0.1),
# # trans = scales::pseudo_log_trans(sigma = 1),
# breaks = c(1, 10, 50, 100, 400), labels = c(1, 10, 50, 100, 400),
# limits = c(0, 500)) + #need to make sure that we include max value of tpm as upper limit
scale_fill_viridis_c(na.value = "transparent") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
strip.text.x = element_text(size=8),
#change the size of the label and the angle so the text is printed horizontally
strip.text.y = element_text(size = 7,angle=0)) -> perma_bbuble
perma_bbuble +
theme(legend.position = "bottom") -> perma_bbuble
perma_bbuble +
scale_y_discrete(label = function(x) stringr::str_trunc(x, 34, side = "center")) ->perma_bbuble
perma_bbuble
perma_bbuble %>%
export::graph2ppt(append = TRUE, width = 20, height = 50,
file = out_pptx)## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxphyseq_AMRgn %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
left_join(tax_mapping,
by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
suffix = c("_x", "")) %>%
mutate(Best_Hit_ARO = fct_reorder(Best_Hit_ARO, Abundance, .desc=FALSE)) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
# df$Item = factor(df$Item, levels = a_lot$Item[order(a_lot$percent)])
# tmp$variable <- factor(melted$variable, levels=rev(levels(melted$variable)))
# mutate(Best_Hit_ARO = fct_reorder(Best_Hit_ARO, logTPM)) %>%
# mutate(Best_Hit_ARO = fct_relevel(Best_Hit_ARO, logTPM)) %>%
# separate_rows(., Drug_Class,sep = '; ',convert = TRUE) %>%
ggplot(mapping = aes(x = as.factor(Day_of_Treatment),
y = Best_Hit_ARO,
fill = logTPM,
)) +
facet_grid(. ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free") + #, labeller = labeller(Drug_Class_multi = label_wrap_gen(width = 15))) +
# facet_grid(. ~ Reactor_Treatment, scales = "free", space = "free", drop = TRUE,labeller = labeller(Drug.Class = label_wrap_gen(width = 15))) +
#we adjust the vertical labels such that it is printed correctly.
geom_tile() +
theme_light() +
# scale_fill_viridis_c(name = "TPM",
# na.value = "transparent", trans = scales::pseudo_log_trans(sigma = 0.1),
# # trans = scales::pseudo_log_trans(sigma = 1),
# breaks = c(1, 10, 50, 100, 400), labels = c(1, 10, 50, 100, 400),
# limits = c(0, 500)) + #need to make sure that we include max value of tpm as upper limit
scale_fill_viridis_c(na.value = "transparent") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
strip.text.x = element_text(size=8),
#change the size of the label and the angle so the text is printed horizontally
strip.text.y = element_text(size = 7,angle=0)) -> perma_bbuble
perma_bbuble +
theme(legend.position = "bottom") +
scale_y_discrete(label = function(x) stringr::str_trunc(x, 34, side = "center")) -> perma_bbuble
perma_bbuble
perma_bbuble %>%
export::graph2ppt(append = TRUE, width = 20, height = 28,
file = out_pptx)## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxphyseq_AMRgn %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
left_join(tax_mapping,
by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
suffix = c("_x", "")) %>%
mutate(Best_Hit_ARO = fct_reorder(Best_Hit_ARO, Abundance, .desc=FALSE)) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
# df$Item = factor(df$Item, levels = a_lot$Item[order(a_lot$percent)])
# tmp$variable <- factor(melted$variable, levels=rev(levels(melted$variable)))
# mutate(Best_Hit_ARO = fct_reorder(Best_Hit_ARO, logTPM)) %>%
# mutate(Best_Hit_ARO = fct_relevel(Best_Hit_ARO, logTPM)) %>%
# separate_rows(., Drug_Class,sep = '; ',convert = TRUE) %>%
ggplot(mapping = aes(x = as.factor(Day_of_Treatment),
y = Best_Hit_ARO,
fill = logTPM,
)) +
facet_grid(Drug_Class_multi ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free") + #, labeller = labeller(Drug_Class_multi = label_wrap_gen(width = 15))) +
# facet_grid(. ~ Reactor_Treatment, scales = "free", space = "free", drop = TRUE,labeller = labeller(Drug.Class = label_wrap_gen(width = 15))) +
#we adjust the vertical labels such that it is printed correctly.
geom_tile() +
theme_light() +
# scale_fill_viridis_c(name = "TPM",
# na.value = "transparent", trans = scales::pseudo_log_trans(sigma = 0.1),
# # trans = scales::pseudo_log_trans(sigma = 1),
# breaks = c(1, 10, 50, 100, 400), labels = c(1, 10, 50, 100, 400),
# limits = c(0, 500)) + #need to make sure that we include max value of tpm as upper limit
scale_fill_viridis_c(na.value = "transparent") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
strip.text.x = element_text(size=8),
#change the size of the label and the angle so the text is printed horizontally
strip.text.y = element_text(size = 7,angle=0)) -> perma_bbuble
perma_bbuble +
theme(legend.position = "bottom") +
scale_y_discrete(label = function(x) stringr::str_trunc(x, 34, side = "center")) -> perma_bbuble
perma_bbuble
perma_bbuble %>%
export::graph2ppt(append = TRUE, width = 20, height = 28,
file = out_pptx)## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptx# perma_bbuble + facet_null() + facet_grid(. ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free") +
# ggside::geom_ysidetile(aes(x = 0, yfill = col[perma_bbuble$data$Drug_Class_multi]), show.legend = TRUE) -> test_annot
#
#
# test_annotpheatmap test:
all
physeq_AMRgn %>%
speedyseq::tax_glom("Best_Hit_ARO") %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
# left_join(tax_mapping,
# by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
# suffix = c("_x", "")) %>%
# mutate(Sample = fct_reorder(Sample, Model,Treatment,Fermentation,Antibiotic_mg.mL, Day_of_Treatment )) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) -> pheatmap_df
pheatmap_df %>%
# replace(is.na(.), 0) %>%
# select(Sample, Best_Hit_ARO, logTPM) %>%
group_by(Sample, Best_Hit_ARO) %>%
pivot_wider(names_from = Sample, Best_Hit_ARO,
values_from = logTPM) %>%
column_to_rownames("Best_Hit_ARO") -> pheatmap_samples
pheatmap_samples %>%
replace(is.na(.), 0) %>%
t() %>%
vegdist(na.rm = TRUE) %>%
hclust() -> clust_t
pheatmap_samples %>%
replace(is.na(.), 0) %>%
vegdist(na.rm = TRUE) %>%
# cor() %>%
hclust() -> clust
pheatmap_samples %>%
pheatmap::pheatmap(cluster_rows = clust ,
cluster_cols = clust_t,
na_col = "transparent",
annotation_row = tax_mapping %>% rownames_to_column('tmp_col') %>%
filter(Best_Hit_ARO %in% taxa_names(physeq_AMRgn)) %>%
distinct(Best_Hit_ARO, .keep_all = TRUE) %>%
column_to_rownames("Best_Hit_ARO") %>%
select(Drug_Class_multi, Resistance_Mechanism_multi),
annotation_colors = list(Drug_Class_multi = col,
Treatment = treat_col),
annotation_col = physeq_AMRgn %>%
sample_data() %>%
data.frame() %>%
select(Model,Treatment,Fermentation,Antibiotic_mg.mL, Day_of_Treatment)) -> p_pheat
p_pheat %>%
print()
p_pheat %>%
export::graph2ppt(append = TRUE, width = 30, height = 30,
file = out_pptx)## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxchick
physeq_AMRgn %>%
speedyseq::tax_glom("Best_Hit_ARO") %>%
subset_samples(Model == "Chicken") %>%
filter_taxa(function(x) sum(x > 0) > 0, TRUE) -> ps_tmp
ps_tmp %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
# left_join(tax_mapping,
# by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
# suffix = c("_x", "")) %>%
# mutate(Sample = fct_reorder(Sample, Model,Treatment,Fermentation,Antibiotic_mg.mL, Day_of_Treatment )) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) -> pheatmap_df
pheatmap_df %>%
# replace(is.na(.), 0) %>%
# select(Sample, Best_Hit_ARO, logTPM) %>%
group_by(Sample, Best_Hit_ARO) %>%
pivot_wider(names_from = Sample, Best_Hit_ARO,
values_from = logTPM) %>%
column_to_rownames("Best_Hit_ARO") -> pheatmap_samples
pheatmap_samples %>%
replace(is.na(.), 0) %>%
t() %>%
vegdist(na.rm = TRUE) %>%
hclust() -> clust_t
# pheatmap_samples %>%
# replace(is.na(.), 0) %>%
# vegdist(na.rm = TRUE) %>%
# # cor() %>%
# hclust() -> clust
pheatmap_samples %>%
pheatmap::pheatmap(cluster_rows = FALSE ,
cluster_cols = clust_t,
na_col = "transparent",
annotation_row = tax_mapping %>% rownames_to_column('tmp_col') %>%
filter(Best_Hit_ARO %in% taxa_names(physeq_AMRgn)) %>%
distinct(Best_Hit_ARO, .keep_all = TRUE) %>%
column_to_rownames("Best_Hit_ARO") %>%
select(Drug_Class_multi, Resistance_Mechanism_multi),
annotation_colors = list(Drug_Class_multi = col)) -> p_pheat #,
# annotation_col = ps_tmp %>%
# sample_data() %>%
# data.frame() %>%
# select(Model,Treatment,Fermentation,Antibiotic_mg.mL, Day_of_Treatment)) -> p_pheat
p_pheat %>%
print()
p_pheat %>%
export::graph2ppt(append = TRUE, width = 20, height = 20,
file = out_pptx)## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxHuman
physeq_AMRgn %>%
speedyseq::tax_glom("Best_Hit_ARO") %>%
subset_samples(Model == "Human") %>%
filter_taxa(function(x) sum(x > 0) > 0, TRUE) -> ps_tmp
ps_tmp %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
# left_join(tax_mapping,
# by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
# suffix = c("_x", "")) %>%
# mutate(Sample = fct_reorder(Sample, Model,Treatment,Fermentation,Antibiotic_mg.mL, Day_of_Treatment )) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) -> pheatmap_df
pheatmap_df %>%
# replace(is.na(.), 0) %>%
# select(Sample, Best_Hit_ARO, logTPM) %>%
group_by(Sample, Best_Hit_ARO) %>%
pivot_wider(names_from = Sample, Best_Hit_ARO,
values_from = logTPM) %>%
column_to_rownames("Best_Hit_ARO") -> pheatmap_samples
pheatmap_samples %>%
replace(is.na(.), 0) %>%
t() %>%
vegdist(na.rm = TRUE) %>%
hclust() -> clust_t
pheatmap_samples %>%
replace(is.na(.), 0) %>%
vegdist(na.rm = TRUE) %>%
# cor() %>%
hclust() -> clust
pheatmap_samples %>%
pheatmap::pheatmap(cluster_rows = clust ,
cluster_cols = clust_t,
na_col = "transparent",
annotation_row = tax_mapping %>% rownames_to_column('tmp_col') %>%
filter(Best_Hit_ARO %in% taxa_names(physeq_AMRgn)) %>%
distinct(Best_Hit_ARO, .keep_all = TRUE) %>%
column_to_rownames("Best_Hit_ARO") %>%
select(Drug_Class_multi, Resistance_Mechanism_multi),
annotation_colors = list(Drug_Class_multi = col,
Treatment = treat_col),
annotation_col = ps_tmp %>%
sample_data() %>%
data.frame() %>%
select(Model,Treatment,Fermentation,Antibiotic_mg.mL, Day_of_Treatment)) -> p_pheat
p_pheat %>%
print()
p_pheat %>%
export::graph2ppt(append = TRUE, width = 20, height = 20,
file = out_pptx)## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptx# see Sneha's & Hannah's code for "plot.df.sep <- separate_rows(plot.df, Drug.Class,sep = '; ',convert = TRUE)"
physeq_AMRgn %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
left_join(tax_mapping,
by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
suffix = c("_x", "")) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
separate_rows(., Resistance_Mechanism,sep = '; ',convert = TRUE) %>%
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Resistance_Mechanism)) +
geom_bar( stat = "identity", colour="black", size=0.025) +
facet_grid(Drug_Class_multi ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("log10 TPM") +
theme_light() -> p_mec
p_mec %>%
ggpubr::set_palette("jco") -> p_mec
p_mec## Warning: Removed 9564 rows containing missing values (position_stack).
p_mec %>%
export::graph2ppt(append = TRUE, width = 20, height = 18,
file = out_pptx)## Warning: Removed 9564 rows containing missing values (position_stack).## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptx# see Sneha's & Hannah's code for "plot.df.sep <- separate_rows(plot.df, Drug.Class,sep = '; ',convert = TRUE)"
physeq_AMRgn %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
left_join(tax_mapping,
by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
suffix = c("_x", "")) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
separate_rows(., Resistance_Mechanism,sep = '; ',convert = TRUE) %>%
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Resistance_Mechanism_multi)) +
geom_bar( stat = "identity", colour="black", size=0.025) +
facet_grid(. ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("log10 TPM") +
theme_light() -> p_mec
p_mec %>%
ggpubr::set_palette("jco") -> p_mec
p_mec## Warning: Removed 9564 rows containing missing values (position_stack).
p_mec %>%
export::graph2ppt(append = TRUE, width = 18, height = 6,
file = out_pptx)## Warning: Removed 9564 rows containing missing values (position_stack).## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptx# see Sneha's & Hannah's code for "plot.df.sep <- separate_rows(plot.df, Drug.Class,sep = '; ',convert = TRUE)"
physeq_AMRgn %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
left_join(tax_mapping,
by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
suffix = c("_x", "")) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
# separate_rows(., Resistance_Mechanism,sep = '; ',convert = TRUE) %>%
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Resistance_Mechanism_multi)) +
geom_bar( stat = "identity", colour="black", size=0.025) +
facet_grid(. ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("log10 TPM") +
theme_light() -> p_mec
p_mec %>%
ggpubr::set_palette("jco") -> p_mec
p_mec## Warning: Removed 9108 rows containing missing values (position_stack).
p_mec %>%
export::graph2ppt(append = TRUE, width = 18, height = 6,
file = out_pptx)## Warning: Removed 9108 rows containing missing values (position_stack).## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptx# see Sneha's & Hannah's code for "plot.df.sep <- separate_rows(plot.df, Drug.Class,sep = '; ',convert = TRUE)"
physeq_AMRgn %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
left_join(tax_mapping,
by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
suffix = c("_x", "")) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
separate_rows(., Resistance_Mechanism,sep = '; ',convert = TRUE) %>%
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Resistance_Mechanism)) +
geom_bar( stat = "identity", position = "fill", colour="black", size=0.025) +
facet_grid(. ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("Proportion") +
theme_light() -> p_mec
p_mec %>%
ggpubr::set_palette("jco") -> p_mec
p_mec## Warning: Removed 9564 rows containing missing values (position_stack).
p_mec %>%
export::graph2ppt(append = TRUE, width = 18, height = 6,
file = out_pptx)## Warning: Removed 9564 rows containing missing values (position_stack).## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptx# https://newbedev.com/how-to-control-ordering-of-stacked-bar-chart-using-identity-on-ggplot2
physeq_AMRgn %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
left_join(tax_mapping,
by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
suffix = c("_x", "")) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
# separate_rows(., Resistance_Mechanism,sep = '; ',convert = TRUE) %>%
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Resistance_Mechanism_multi)) +
geom_bar( stat = "identity", position = "fill", colour="black", size=0.00125) +
facet_grid(. ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("Proportion") +
theme_light() -> p_mec
p_mec %>%
ggpubr::set_palette("jco") -> p_mec
p_mec## Warning: Removed 9108 rows containing missing values (position_stack).
p_mec %>%
export::graph2ppt(append = TRUE, width = 18, height = 6,
file = out_pptx)## Warning: Removed 9108 rows containing missing values (position_stack).## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxphyseq_drug_multi %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Drug_Class_multi)) +
geom_bar( stat = "identity", colour="black", size=0.025) +
scale_fill_manual(values=col) +
#geom_bar(stat="identity",colour=NA,size=0) +
facet_grid( ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("log10 TPM") +
theme_light() -> p_drug
p_drug
p_drug %>%
export::graph2ppt(append = TRUE, width = 18, height = 6,
file = out_pptx)## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxphyseq_drug_multi %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Drug_Class_multi)) +
geom_bar( stat = "identity", position="fill", colour="black", size=0.025) +
scale_fill_manual(values=col) +
#geom_bar(stat="identity",colour=NA,size=0) +
facet_grid(. ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x", drop = TRUE) +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("Proportion") +
theme_light() -> p_drug_prop
p_drug_prop
p_drug_prop %>%
export::graph2ppt(append = TRUE, width = 18, height = 6,
file = out_pptx)## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxphyseq_AMRgn_round <- physeq_AMRgn
round(otu_table(physeq_AMRgn_round)) -> otu_table(physeq_AMRgn_round)
physeq_AMRgn_round %>%
phyloseq_alphas(phylo = FALSE) -> alpha_df## Loading required package: metagMisc##
## Attaching package: 'metagMisc'## The following object is masked from 'package:purrr':
##
## some## Loading required package: microbiome##
## microbiome R package (microbiome.github.com)
##
##
##
## Copyright (C) 2011-2020 Leo Lahti,
## Sudarshan Shetty et al. <microbiome.github.io>##
## Attaching package: 'microbiome'## The following object is masked from 'package:scales':
##
## alpha## The following object is masked from 'package:vegan':
##
## diversity## The following object is masked from 'package:ggplot2':
##
## alpha## The following object is masked from 'package:base':
##
## transform## Observed richness## Other forms of richness## Diversity## Evenness## Dominance## Rarity## New names:
## * id -> id...1
## * id -> id...62
## * id -> id...85# plot panels for each treatment
alpha_df %>%
pivot_longer(cols = all_of(c("Observed", "diversity_shannon")), values_to = "value", names_to = 'alphadiversiy', values_drop_na = TRUE) %>%
mutate(alphadiversiy = fct_relevel(alphadiversiy, c("Observed", "diversity_shannon"))) %>%
# filter(!is.na(Fermentation)) %>%
# filter(Treatment != "DONOR") %>%
ggplot(aes(x = Day_of_Treatment, y = value, color = Treatment)) +
geom_point() +
geom_line(linetype = "dashed", size = 0.25) +
labs(title = "AMR Gene Richness",
y = "Richness", x = "Days of Treatment") +
scale_color_manual(values = treat_col,
na.value = "black") +
facet_grid(alphadiversiy ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "fixed", drop = TRUE) +
theme_light() -> p_alpha
p_alpha## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
p_alpha %>%
export::graph2ppt(append = TRUE, width = 20, height = 6,
file = out_pptx)## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxBeta Diversity:
# sample_data(physeq)$Reactor_Treatment <- fct_relevel(sample_data(physeq)$Reactor_Treatment, "DONOR", "CR_UNTREATED", "TR1_CTX+HV292.1",
# "TR2_CTX","TR3_HV292.1", "TR4_VAN", "TR5_VAN+CCUG59168", "TR6_CCUG59168")
#
# sample_data(physeq)$Treatment <- fct_relevel(sample_data(physeq)$Treatment, "DONOR", "UNTREATED", "CTX+HV292.1", "CTX","HV292.1","VAN+CCUG59168", "VAN", "CCUG59168")
# physeq %>%
# rarefy_even_depth(sample.size = 43,
# rngseed = 123) -> phyloseq_rare
physeq_AMRgn %>%
phyloseq_compute_bdiv(norm = "pc",
phylo = FALSE,
seed = 123) -> bdiv_list## Loading required package: ape## Loading required package: GUniFrac## Registered S3 method overwritten by 'rmutil':
## method from
## print.response httrphyseq_AMRgn %>%
# subset_samples(Treatment != "DONOR") %>%
phyloseq_plot_bdiv(dlist = bdiv_list,
seed = 123,
axis1 = 1,
axis2 = 2) -> pcoa## [1] "bray"
## [1] "sorensen"
## [1] "bjaccard"
## [1] "wjaccard"# phyloseq_plot_bdiv(bdiv_list,
# # m = "CoDa",
# seed = 123) -> coda
#
pcoa$wjaccard$layers = NULL
pcoa$wjaccard + geom_point(size=3,
aes(color = Treatment,
fill = NULL,
shape = Fermentation %>% as.factor(),
alpha = Day_of_Treatment)) +
theme_light() +
geom_path(arrow = arrow(type = "open", angle = 30, length = unit(0.15, "inches")),
size = 0.05, linetype = "dashed", inherit.aes = TRUE, aes(group=Reactor_Treatment), show.legend = FALSE) +
scale_alpha_continuous(range=c( 0.9, 0.3)) +
scale_fill_manual(values = treat_col,
na.value = "black") +
scale_color_manual(values = treat_col,
na.value = "black") +
scale_shape_manual(values = c(15, 19), na.value = 17) +
theme_classic() -> p1 # +
# labs(col=NULL, fill = NULL, shape = NULL) + guides(shape=TRUE) -> p1
p1 + facet_grid(Model ~.)
p1 %>%
export::graph2ppt(append = TRUE, width = 8, height = 6,
file = out_pptx)## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxphyseq_AMRgn %>%
tax_table() %>%
data.frame() %>%
rownames_to_column("id") %>%
mutate(gene = id) %>%
column_to_rownames("id") %>%
as.matrix %>%
tax_table() -> tax_table(physeq_AMRgn)
physeq_AMRgn %>%
# subset_samples(Treatment != "DONOR") %>%
phyloseq_add_taxa_vector(dist = bdiv_list$wjaccard,
phyloseq = .,
figure_ord = p1,
tax_rank_plot = "Best_Hit_ARO", taxrank_glom = "Best_Hit_ARO",
top_r = 10, fact = 0.6) -> pco_env
# pco_env$plot %>%
# export::graph2ppt(append = TRUE, width = 8, height = 6,
# file = out_pptx)
pco_env$signenvfit %>%
DT::datatable()physeq_drug_multi %>%
# subset_samples(Treatment != "DONOR") %>%
phyloseq_add_taxa_vector(dist = bdiv_list$wjaccard,
phyloseq = .,
figure_ord = p1,
tax_rank_plot = "Drug_Class_multi", taxrank_glom = "Drug_Class_multi",
top_r = 10, fact = 0.6) -> pco_env
pco_env$plot + facet_grid(Model ~ .) -> pco_env$plot
pco_env$plot## Warning: Removed 2 rows containing missing values (geom_text_repel).
pco_env$plot %>%
export::graph2ppt(append = TRUE, width = 8, height = 6,
file = out_pptx)## Warning: Removed 2 rows containing missing values (geom_text_repel).## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxpco_env$signenvfit %>%
DT::datatable()Make the heatmap pretty
physeq %>%
# speedyseq::tax_glom("Best_Hit_ARO") %>%
speedyseq::psmelt() %>%
rename(TPM = Abundance) %>%
mutate(TPM = na_if(TPM, 0)) %>%
# mutate(logTPM = log10(TPM + 1)) %>%
# filter(Drug_Class %in% c("cephalosporin", "glycopeptide")) %>%
# mutate(Drug_Class = fct_relevel(Drug_Class, c("cephalosporin", "glycopeptide"))) %>%
filter(Best_Hit_ARO %in% c("vanZA", "vanYA", "vanXA", "vanA", "vanHA", "vanRA", "CTX-M-1")) %>%
# mutate(Classification = fct_relevel(Classification, c("MAG", "Chromosome", "Plasmid", "Uncertain - plasmid or chromosomal"))) %>%
ggplot(mapping = aes(x = Day_of_Treatment ,
y = Best_Hit_ARO ,
fill = TPM,
)) +
facet_grid(Drug_Class ~ Model + Treatment + Fermentation + Antibiotic_mg.mL, scales = "free", space = "fixed", drop = TRUE) +
# facet_grid(. ~ Treatment, scales = "free", space = "free", drop = TRUE) +
# facet_grid(. ~ Reactor_Treatment, scales = "free", space = "free", drop = TRUE,labeller = labeller(Drug_Class = label_wrap_gen(width = 15))) +
#we adjust the vertical labels such that it is printed correctly.
geom_tile() +
theme_bw() +
# scale_fill_viridis_c(name = "TPM",
# na.value = "transparent", trans = scales::pseudo_log_trans(sigma = 0.1),
# # trans = scales::pseudo_log_trans(sigma = 1),
# breaks = c(1, 10, 50, 100, 400), labels = c(1, 10, 50, 100, 400),
# limits = c(0, 500)) + #need to make sure that we include max value of tpm as upper limit
scale_fill_viridis_c(na.value = "transparent") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
strip.text.x = element_text(size=8),
#change the size of the label and the angle so the text is printed horizontally
strip.text.y = element_text(size = 7,angle=0)) -> perma_bbuble
perma_bbuble +
theme(legend.position = "bottom") -> perma_bbuble
perma_bbuble
perma_bbuble %>%
export::graph2ppt(append = TRUE, width = 16, height = 4,
file = out_pptx)## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxWhat’s next? - MAG, plasmid, Chromosome … proportions. - SQM contig taxonomy coloring/facetting.
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_normalisation.R")
ps_tmp_tax <- physeq
ps_tmp_tax %>%
tax_table() %>%
data.frame() %>%
# select(Tax) %>%
separate(Tax, into = c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"), sep = ";") %>%
mutate(across(c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"), gsub, pattern = "[a-z]_", replacement = "")) %>%
as.matrix() -> tax_table(ps_tmp_tax)## Warning: Expected 7 pieces. Additional pieces discarded in 23 rows [8, 10,
## 11, 13, 41, 44, 58, 79, 83, 95, 110, 135, 138, 139, 166, 181, 189, 207, 210,
## 220, ...].## Warning: Expected 7 pieces. Missing pieces filled with `NA` in 111 rows [4, 6,
## 9, 16, 17, 21, 23, 26, 29, 30, 33, 36, 45, 48, 51, 52, 53, 55, 56, 57, ...].ps_tmp_tax_2 <- ps_tmp_tax
ps_tmp_tax %>%
tax_table() %>%
data.frame() %>%
select(c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")) %>%
as.matrix() -> tax_table(ps_tmp_tax_2)
ps_tmp_tax_2 %>%
phyloseq_get_strains_fast(species = FALSE) %>%
tax_table() %>%
data.frame() %>%
mutate(Strain = gsub(" human_megahit.*","",Strain)) %>%
mutate(Strain = gsub(" chick_megahit.*","",Strain)) -> full_tax## Joining, by = "ASV"ps_tmp_tax %>%
tax_table() %>%
data.frame() %>%
rownames_to_column("ASV") %>%
# select(-Tax) %>%
left_join(full_tax %>%
select(Strain) %>% rownames_to_column("ASV"),
by = c("ASV" = "ASV")) %>%
select(-Predicted_DNA) %>%
column_to_rownames("ASV") %>%
as.matrix() -> tax_table(ps_tmp_tax)
tax_table(ps_tmp_tax) %>%
data.frame() %>%
select(c("Best_Hit_ARO", "Strain")) %>%
head()## Best_Hit_ARO Strain
## human_megahit_1512_86-616 emrR Escherichia coli
## human_megahit_1512_743-1915 emrA Escherichia coli
## human_megahit_1512_1932-3470 emrB Escherichia coli
## human_megahit_1935_13232-13930 vanRD unknown Firmicutes (Class)
## human_megahit_4077_533-1408 CTX-M-1 Escherichia coli
## human_megahit_4631_5392-6930 emrB unknown Citrobacter (Genus)ps_tmp_tax %>%
saveRDS("~/Projects/ETH/Alessia/mobilome/ps_combined_tax_clean.rds")
as(tax_table(ps_tmp_tax), "matrix") %>%
as.data.frame() %>%
rownames_to_column('ASV') -> tax_table_ps_combined
as(otu_table(ps_tmp_tax), "matrix") %>%
as.data.frame() %>%
rownames_to_column('ASV') -> otu_table_ps_combined
ps_tmp_tax %>%
sample_data() %>%
data.frame() %>%
rownames_to_column('sample_id') -> meta_ps_combined
otu_table_ps_combined %>%
# left_join(refseq_df, by = 'ASV') %>%
left_join(tax_table_ps_combined, by = c("ASV" = "ASV")) %>%
dplyr::select(ASV, everything()) -> merged_ps_combined
merged_ps_combined %>%
write_tsv("~/Projects/ETH/Alessia/mobilome/ps_combined_tax_clean_rds.tsv")tax_table(ps_tmp_tax) %>%
data.frame() %>%
select(c("Best_Hit_ARO", "Strain")) %>%
group_by(Strain) %>%
add_count() %>%
arrange(-n)## # A tibble: 244 × 3
## # Groups: Strain [29]
## Best_Hit_ARO Strain n
## <chr> <chr> <int>
## 1 emrR Escherichia coli 99
## 2 emrA Escherichia coli 99
## 3 emrB Escherichia coli 99
## 4 CTX-M-1 Escherichia coli 99
## 5 baeS Escherichia coli 99
## 6 CRP Escherichia coli 99
## 7 YojI Escherichia coli 99
## 8 evgS Escherichia coli 99
## 9 evgA Escherichia coli 99
## 10 Escherichia coli mdfA Escherichia coli 99
## # … with 234 more rowsDefine colors for Taxonomic information:
ps_tmp_tax %>% speedyseq::psmelt() %>% pull(Strain) %>% unique() -> strains
set.seed(456)
col_strains <- distinctColorPalette(length(strains),
runTsne = FALSE)
names(col_strains) <- strains# see Sneha's & Hannah's code for "plot.df.sep <- separate_rows(plot.df, Drug.Class,sep = '; ',convert = TRUE)"
ps_tmp_tax %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
# left_join(tax_mapping,
# by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
# suffix = c("_x", "")) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
# separate_rows(., Resistance_Mechanism,sep = '; ',convert = TRUE) %>%
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Strain)) +
geom_bar( stat = "identity", colour="black", size=0.025) +
facet_grid(Drug_Class_multi ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("log10 TPM") +
theme_light() +
scale_fill_manual(name = "", values = col_strains) +
theme(legend.position = "bottom") -> p_mec
# p_mec %>%
# ggpubr::set_palette("jco") -> p_mec
p_mec## Warning: Removed 15421 rows containing missing values (position_stack).
p_mec %>%
export::graph2ppt(append = TRUE, width = 20, height = 18,
file = out_pptx)## Warning: Removed 15421 rows containing missing values (position_stack).## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxps_tmp_tax %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
# left_join(tax_mapping,
# by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
# suffix = c("_x", "")) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
# separate_rows(., Resistance_Mechanism,sep = '; ',convert = TRUE) %>%
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Drug_Class_multi)) +
geom_bar( stat = "identity", colour="black", size=0.025) +
facet_grid( Strain ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("log10 TPM") +
theme_light() +
scale_fill_manual(name = "", values = col) +
theme(legend.position = "bottom") -> p_mec
# p_mec %>%
# ggpubr::set_palette("jco") -> p_mec
p_mec## Warning: Removed 15421 rows containing missing values (position_stack).
p_mec %>%
export::graph2ppt(append = TRUE, width = 20, height = 18,
file = out_pptx)## Warning: Removed 15421 rows containing missing values (position_stack).## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptx# see Sneha's & Hannah's code for "plot.df.sep <- separate_rows(plot.df, Drug.Class,sep = '; ',convert = TRUE)"
ps_tmp_tax %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
# left_join(tax_mapping,
# by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
# suffix = c("_x", "")) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
# separate_rows(., Resistance_Mechanism,sep = '; ',convert = TRUE) %>%
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Strain)) +
geom_bar( stat = "identity", colour="black", position = "fill", size=0.05) +
facet_grid(Drug_Class_multi ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("Proportion") +
theme_light() +
scale_fill_manual(name = "", values = col_strains) +
theme(legend.position = "bottom") -> p_mec
# p_mec %>%
# ggpubr::set_palette("jco") -> p_mec
p_mec## Warning: Removed 15421 rows containing missing values (position_stack).
p_mec %>%
export::graph2ppt(append = TRUE, width = 20, height = 18,
file = out_pptx)## Warning: Removed 15421 rows containing missing values (position_stack).## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptx# see Sneha's & Hannah's code for "plot.df.sep <- separate_rows(plot.df, Drug.Class,sep = '; ',convert = TRUE)"
ps_tmp_tax %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
# left_join(tax_mapping,
# by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
# suffix = c("_x", "")) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
# separate_rows(., Resistance_Mechanism,sep = '; ',convert = TRUE) %>%
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Strain)) +
geom_bar( stat = "identity", colour="black", size=0.025) +
facet_grid(. ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("log10 TPM") +
theme_light() +
scale_fill_manual(name = "", values = col_strains) +
theme(legend.position = "bottom") -> p_mec
# p_mec %>%
# ggpubr::set_palette("jco") -> p_mec
p_mec## Warning: Removed 15421 rows containing missing values (position_stack).
p_mec %>%
export::graph2ppt(append = TRUE, width = 18, height = 6,
file = out_pptx)## Warning: Removed 15421 rows containing missing values (position_stack).## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptx# see Sneha's & Hannah's code for "plot.df.sep <- separate_rows(plot.df, Drug.Class,sep = '; ',convert = TRUE)"
ps_tmp_tax %>%
microbiome::transform(transform = "log10p") %>%
speedyseq::psmelt() %>%
# left_join(tax_mapping,
# by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
# suffix = c("_x", "")) %>%
mutate(Abundance = na_if(Abundance, 0)) %>%
mutate(logTPM = log10(Abundance + 1)) %>%
# separate_rows(., Resistance_Mechanism,sep = '; ',convert = TRUE) %>%
# ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = fct_reorder(Strain, logTPM))) +
ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = reorder(Strain, logTPM, FUN=median, na.rm = TRUE))) +
geom_bar( stat = "identity", position="fill", colour="black", size=0.025) +
facet_grid(. ~ Model + Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
ggtitle("ARG Type Abundance per Sample") +
xlab("Day (Treatment)") +
ylab("Proportion") +
theme_light() +
scale_fill_manual(name = "", values = col_strains) +
theme(legend.position = "bottom") -> p_mec
# p_mec %>%
# ggpubr::set_palette("jco") -> p_mec
p_mec## Warning: Removed 15421 rows containing missing values (position_stack).
p_mec %>%
export::graph2ppt(append = TRUE, width = 18, height = 6,
file = out_pptx)## Warning: Removed 15421 rows containing missing values (position_stack).## Exported graph as ~/Projects/ETH/Alessia/mobilome/graph.pptxsessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] GUniFrac_1.4 ape_5.6 microbiome_1.14.0
## [4] metagMisc_0.0.4 gdtools_0.2.3 reshape2_1.4.4
## [7] scales_1.1.1 randomcoloR_1.1.0.1 vegan_2.5-7
## [10] lattice_0.20-45 permute_0.9-5 RColorBrewer_1.1-2
## [13] microViz_0.9.0 here_1.0.1 ggrepel_0.9.1
## [16] speedyseq_0.5.3.9018 phyloseq_1.36.0 forcats_0.5.1
## [19] stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4
## [22] readr_2.1.0 tidyr_1.1.4 tibble_3.1.6
## [25] ggplot2_3.3.5 tidyverse_1.3.1.9000
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.1 lme4_1.1-27.1
## [4] htmlwidgets_1.5.4 grid_4.1.2 Rtsne_0.15
## [7] modeest_2.4.0 munsell_0.5.0 codetools_0.2-18
## [10] statmod_1.4.36 DT_0.20 withr_2.4.3
## [13] colorspace_2.0-2 Biobase_2.52.0 highr_0.9
## [16] knitr_1.36 uuid_1.0-3 rstudioapi_0.13
## [19] stats4_4.1.2 ggsignif_0.6.3 officer_0.4.0
## [22] labeling_0.4.2 GenomeInfoDbData_1.2.6 bit64_4.0.5
## [25] farver_2.1.0 pheatmap_1.0.12 rhdf5_2.36.0
## [28] fBasics_3042.89.1 rprojroot_2.0.2 vctrs_0.3.8
## [31] generics_0.1.1 xfun_0.28 R6_2.5.1
## [34] GenomeInfoDb_1.28.4 clue_0.3-60 bitops_1.0-7
## [37] rhdf5filters_1.4.0 assertthat_0.2.1 vroom_1.5.6
## [40] gtable_0.3.0 spatial_7.3-14 timeDate_3043.102
## [43] rlang_0.4.12 systemfonts_1.0.2 splines_4.1.2
## [46] rstatix_0.7.0 broom_0.7.11 rgl_0.107.14
## [49] yaml_2.2.1 abind_1.4-5 modelr_0.1.8
## [52] crosstalk_1.2.0 backports_1.4.1 tools_4.1.2
## [55] ellipsis_0.3.2 jquerylib_0.1.4 biomformat_1.20.0
## [58] BiocGenerics_0.38.0 stargazer_5.2.2 stabledist_0.7-1
## [61] devEMF_4.0-2 Rcpp_1.0.7 plyr_1.8.6
## [64] base64enc_0.1-3 zlibbioc_1.38.0 RCurl_1.98-1.5
## [67] ggpubr_0.4.0 rpart_4.1-15 viridis_0.6.2
## [70] statip_0.2.3 S4Vectors_0.30.2 haven_2.4.3
## [73] cluster_2.1.2 fs_1.5.2 magrittr_2.0.1
## [76] data.table_1.14.2 timeSeries_3062.100 openxlsx_4.2.4
## [79] lmerTest_3.1-3 flextable_0.6.9 reprex_2.0.1
## [82] matrixStats_0.61.0 hms_1.1.1 evaluate_0.14
## [85] xtable_1.8-4 rio_0.5.27 readxl_1.3.1
## [88] IRanges_2.26.0 gridExtra_2.3 compiler_4.1.2
## [91] V8_3.6.0 crayon_1.4.2 minqa_1.2.4
## [94] htmltools_0.5.2 mgcv_1.8-38 tzdb_0.2.0
## [97] lubridate_1.8.0 DBI_1.1.1 rmutil_1.1.5
## [100] dbplyr_2.1.1 MASS_7.3-54 boot_1.3-28
## [103] Matrix_1.3-4 ade4_1.7-18 car_3.0-11
## [106] cli_3.1.0 parallel_4.1.2 igraph_1.2.10
## [109] pkgconfig_2.0.3 numDeriv_2016.8-1.1 foreign_0.8-81
## [112] xml2_1.3.2 foreach_1.5.1 bslib_0.3.1
## [115] multtest_2.48.0 XVector_0.32.0 rvg_0.2.5
## [118] rvest_1.0.2 digest_0.6.29 Biostrings_2.60.2
## [121] rmarkdown_2.11 cellranger_1.1.0 curl_4.3.2
## [124] export_0.3.0 nloptr_1.2.2.2 lifecycle_1.0.1
## [127] nlme_3.1-153 jsonlite_1.7.2 Rhdf5lib_1.14.2
## [130] carData_3.0-4 viridisLite_0.4.0 fansi_0.5.0
## [133] pillar_1.6.4 ggsci_2.9 fastmap_1.1.0
## [136] httr_1.4.2 survival_3.2-13 glue_1.6.0
## [139] zip_2.2.0 iterators_1.0.13 bit_4.0.4
## [142] stringi_1.7.6 sass_0.4.0 stable_1.1.4